v20200426.2342

Objetivos

  • Identificar e sumariar dados categóricos em tabela de contingência.
  • Identificar distribuição amostral \(\chi^2\) (qui-quadrado).
  • Formular hipótese estatística na análise de dados categóricos.
  • Definir e determinar frequências observada e esperada.
  • Identificar, aplicar e interpretar os testes de aderência.
  • Identificar e formular análise de variáveis categóricas em tabela de contingência \(l\)x\(c\) (números de linhas versus de colunas).
  • Analisar tabelas de contingencia com variáveis ordinais.
  • Analisar de tabelas “tridimensionais” \(l\)x\(c\)x\(k\), com a inclusão de variável de estratificação.
  • Realizar os procedimentos estatísticos em R.

Roteiro

Há diversas funções e procedimentos envolvidos neste texto. Como roteiro orientador, verifique:

Delineamento entre participantes

teste de independência / dissociação

  • Tabela unidimensional
    • Uma variável nominal:
      • Teste de aderência: chisq.test
  • Tabela bidimensional
    • Duas variáveis nominais dicotômicas
      • Teste qui-quadrado de Pearson: chisq.test
      • Teste exato de Fisher de OR bruto: fisher.test
      • Teste qui-quadrado mínimo
    • Uma variável nominal dicotômica e uma ordinal
      • Teste qui-quadrado de Cochran-Armitage (trend): - - DescTools::CochranArmitageTest
    • Uma variável nominal politômica e uma ordinal
      • Teste qui-quadrado de Cochran-Armitage generalizado: coin::chisq_test
    • Duas variáveis ordinais
      • Teste qui-quadrado de Mantel-Haenszel: DescTools::MHChisqTest
  • Tabela tridimensional
    • Duas variáveis nominais dicotômicas e uma variável de estratificação nominal politômica
      • Teste de Mantel-Haenszel de OR comum: stats::mantelhaen.test
    • Duas variáveis nominais politômicas e uma variável de estratificação nominal politômica
      • Teste de Mantel-Haenszel de OR comum generalizado: coin::cmh_test
    • Duas variáveis ordinais e uma variável de estratificação nominal politômica
      • Teste qui-quadrado ordinal estratificado (linear-by-linear association model): coin::lbl_test

Delineamento intraparticipantes

teste de concordância entre dois métodos de avaliação

  • Tabela bidimensional
    • Duas variáveis nominais dicotômicas
      • Teste de McNemar: mcnemar.test, coin::mh_test
      • Teste kappa de Cohen: epiR::epi.kappa

Contingência

Em sentido geral, contingência pode significar qualquer relação de dependência entre eventos ambientais ou entre eventos comportamentais e ambientais (Catania, 1993; Skinner, 1953, 1969; Todorov, 1985). Embora possa ser encontrado nos dicionários com diferentes significados, esse termo é empregado, na análise do comportamento, como termo técnico para enfatizar como a probabilidade de um evento pode ser afetada ou causada por outros eventos (Catania, 1993, p. 368).
Deisy das Graças de Souza (UFSCar)
http://www.itcrcampinas.com.br/txt/texto_deisy.pdf

\[~\]

Relações entre duas variáveis

1
Exposição e Desfecho

tabela 2x2

independência

concordância

discordância

\[~\]

Associação

2
Não avalia
causa e efeito:
apenas associação.

causa e efeito positivo?

causa e efeito negativo?

causa e efeito positivo?

Exemplos

Gatos e alergia

- estresse

- mecanismo imune

- desinfetantes

- alergia e estresse

- alergia e gatos

Vinho e doença coronariana

- ambiente inseguro

- moradia precária

- má nutrição

- boa moradia

- segurança

- boa alimentação

- outros cuidados

- é o vinho?

Fumo e Câncer de Pulmão

- The Cigarette Century

- How could one attribute lung cancer to cigarette smoking?

“The public impression of scientific and medical uncertainty that resulted was a crucial element in maintaining the market of current smokers as well as recruiting new ones. Industry literature, for example, frequently pointed to the fact that nonsmokers could also develop lung cancer. Therefore, they argued, how could one attribute lung cancer to cigarette smoking?

- Medical knowledge is provisional

medical knowledge was always provisional and contingent. Just as drugs deemed “effective” do not work in every case, so too a cause of disease does not always result in disease. As Richard Doll would later explain, the epidemiologists had identified a cause of lung cancer (and other diseases), not the cause.

- Fisher was against the association between smoking and lung cancer

Another vocal critic of the lung cancer findings was Sir Ronald Fisher, the leading biometrician and geneticist in Great Britain during the first half of the twentieth century and a man deeply committed to bringing statistical analysis to genetics and agricultural experimentation. His 1925 book, Statistical Methods for Research Workers, quickly became a classic, leading to academic appointments at University College London and Cambridge University. Fisher’s critiques were similar to Berkson’s. The ethical impossibility of conducting a randomized experiment led him to question the results of the epidemiological studies. As a believer in genetic* notions of cancer causality, Fisher speculated that there was some constitutional factor that led individuals both to become smokers and to get lung cancer, even though smoking and lung cancer might not be causally related. Doll and Hill repeatedly rebutted this theory, returning to the critical question of how to account for the rise in lung cancers during the twentieth century if the disease was simply “constitutional.”

* 1923: One of the [A. Bradford] Hill’s innovations was the first randomized, double-blind clinical trial, designed to reduce investigator bias in the evaluation of clinical outcomes. […] This method, which drew on Fisher’s agricultural experimentation in genetics, became a critical new tool for evaluating medical treatments.

By the 1930s, the relationship of smoking to cancer was a topic of uresolved debate. It was the life insurance industry, which like the tobacco corporations grew by leaps and bounds in the first half of the twentieth century, that took the lead in understanding the effects of smoking on health.

By the early 1950s, however, it was abundantly clear that the evidence implicating cigarette smoking as a risk to health was now of a different order. First, the link between smoking and disease was categorical, outside the realm of individual clinical judgment.

They [tobacco industry] responded with a new and unprecedent public relations strategy. Its goals was to produce and sustain scientific skepticism and controversy in order to disrupt the emerging consensus on the harms of cigarette smoking. This strategy required intrusions into scientific process and procedure. […] So long as there appeared to be doubt, so long as the industry could assert “not proven”. […] The future of the cigarette would now depend on the successful production of a scientific controversy.

Óculos e xixi

Quando alguém lhe traz uma opinião coincidente com sua crença, você a toma como um fato. Quando alguém lhe traz um fato discordante de sua crença, você o toma como uma opinião.

Prof. Eduardo Massad

A hipótese nula

3
Tomada de decisão:
Hipótese nula (\(H_0\))
e hipótese alternativa (\(H_1\))

Adotaremos, neste texto, a seguinte convenção:

Desfecho + Desfecho -
Exposição + a b
Exposição - c d

Teste Qui-quadrado

Os testes baseados na estatística qui-quadrado podem ser divididos em dois tipos:

  1. Teste de aderência uma variável nominal com duas ou mais categorias: teste de frequências hipotetizadas

  2. Teste de independência entre duas variáveis nominais com duas ou mais categorias

- a distribuição \(\chi^2\)

Antes de discutirmos os testes, vamos explorar um pouco a distribuição que utilizam. Diferente da distribuição normal que têm domínio de \([- \infty, + \infty]\), estas são distribuições assimétricas, iniciadas em zero, que tendem assintoticamente a zero com o valor crescente de \(\chi^2\), \([0, + \infty]\).

O comando para calcular a probabilidade acumulada entre zero e cada valor de qui-quadrado é:

dchisq(x, df, ncp = 0, log = FALSE)

onde x é o valor de \(\chi^2\), df são os graus de liberdade (degrees of freedom) e ncp é o parâmetro de não centralidade (non-centrality parameter = 0, por default, correspondente à distribuição centrada de \(H_0\); voltaremos ao ncp adiante). O parâmetro log, desligado por default é para a possibilidade de devolver o logarítmo do valor (não o exploraremos neste texto).

Pode-se gerar curvas desta distribuição com o seguinte código, mostrando seu aspecto entre 2 e 6 graus de liberdade:

chi2.valor <- seq(from=0, to=40, by=0.01)
graus.liberdade <- 2:6
for (g in 1:length(graus.liberdade))
{
  p <- dchisq(chi2.valor, df=graus.liberdade[g])
  plot(chi2.valor, p, type="l", 
       main=paste("Distribuição de Qui^2, df=",
                  graus.liberdade[g],sep=""),
       xlab="Qui-quadrado", ylab="Densidade")
}

ou, caso prefira ver todas as curvas no mesmo gráfico, sugere-se:

source("friendlycolor.R")
cor <- 3
padrao <- 1
leg.txt <- c()
leg.cor <- c()
chi2.valor <- seq(from=0, to=40, by=0.01)
graus.liberdade <- 2:6
for (g in 1:length(graus.liberdade))
{
  p <- dchisq(chi2.valor, df=graus.liberdade[g])
  if (g==1)
  {
    plot(chi2.valor, p, type="l", lty=padrao,
         col=friendlycolor(cor), lwd=2,
         main="Distribuições de Qui^2",
         xlab="Qui-quadrado", ylab="Densidade")
    
  } else
  {
    lines(chi2.valor, p, lty=padrao,
         col=friendlycolor(cor), lwd=2)
  }
  # guarda para a legenda
  leg.txt <- c(leg.txt, paste("df = ",graus.liberdade[g],sep=""))
  leg.cor <- c(leg.cor, friendlycolor(cor))
  cor <- cor+6
  padrao <- padrao+1
}
legend("topright",
       leg.txt, 
       col=leg.cor,
       lwd=2, 
       lty=1:5, 
       box.lwd=0, bg="transparent")  

Teste de aderência: teste qui-quadrado de uma variável

Permite descobrir se um conjunto de frequências observadas difere de um outro conjunto de frequências hipotetizadas.

suposição

Independência das observações: cada observação da unidade observacional deve ser alocada em apenas um categoria da variável nominal.
SIEGEL, S. & CASTELLAN Jr., N.J. (1988) Nonparametric statistics for the behavioral sciences. 2ª ed. NY: McGraw-Hill, p. 111

exemplo 1: chocolate

Cento e dez pessoas foram solicitadas a manifestar suas preferências com respeito a chocolates. Queremos verificar se alguma das marcas é preferida em detrimento de outras.

Se não são (\(H_0\)), devemos esperar aproximadamente o mesmo número de pessoas em cada categoria.

Com certeza, não haverá exatamente o mesmo número de participantes em cada categoria, mas os números devem ter uma distribuição com razoável proximidade.

O número de pessoas que escolheu apenas uma entre quatro marcas diferentes está em chocolate.xlsx.

O teste está implementado em TesteQuiQuadradoAderencia_Chocolate.R:

# TesteQuiQuadradoAderencia_Chocolate.R

library(readxl)

cat("Teste de aderencia a distribuicao multinomial\n")

df_chocolate <- read_excel("chocolate.xlsx")

prmatrix(df_chocolate, rowlab = rep("",nrow(df_chocolate)))

out <- chisq.test(df_chocolate$Pessoas, p=df_chocolate$probH0, 
                  simulate.p.value = TRUE, B = 1e6) 
df <- chisq.test(df_chocolate$Pessoas, p=df_chocolate$probH0)$parameter
X2 <- out$statistic
p <- out$p.value
n <- sum(chisq.test(df_chocolate$Pessoas, p=df_chocolate$probH0)$observed)
V <- sqrt((X2/n)/df)
cat("X^2 =", X2 , "gl =",  df, "p exato =", p, "\n")
cat("\nHeuristica de significancia: X^2/gl > 2:\n", X2/df, "\n")
cat("\nResiduos estandardizados:\n", out$stdres, "\n")
if (0 <= V & V < 0.1) {gV <- "minimo"}
if (0.1 <= V & V < 0.3) {gV <- "pequeno"}
if (0.3 <= V & V < 0.5) {gV <- "intermediario"}
if (0.5 <= V & V <= 1.0) {gV <- "grande"}
cat("V de Cramer = ", V,"\n",sep="")
cat("(Grau ", gV, "de dessemelhanca entre as distribuicoes observada e hipotetizada\n",sep="")

Observe o código:

  • lê a planilha no formato xlsx com os dados e a exibe.
  • note que denominamos uma coluna como probH0, contendo \(0.25\) em todas as linhas; estas são as frequências esperadas para a hipótese nula, \(H_0\).
  • O teste em si recebe a coluna com o número de pessoas que preferiram cada chocolate e estas frequências esperadas; é executado com bootstrapping, com \(10^5\) iterações (parâmetro \(B\)).
  • O restante do código serve para reorganizar e ecoar os resultados, incluindo o cálculo de tamanho de efeito dado pelo V de Cramér.

Este código gera a seguinte saída:

Teste de aderencia a distribuicao multinomial
 Chocolate Pessoas probH0
 "A"       "20"    "0.25"
 "B"       "60"    "0.25"
 "C"       "10"    "0.25"
 "D"       "20"    "0.25"
X^2 = 53.63636 gl = 3 p exato = 9.99999e-07 

Heuristica de significancia: X^2/gl > 2:
 17.87879 

Residuos estandardizados:
 -1.651446 7.156264 -3.853373 -1.651446 
V de Cramer = 0.4031556
(Grau intermediariode dessemelhanca entre as distribuicoes observada e hipotetizada

Com \(p=9.999 \cdot 10^-5\), rejeita-se \(H_0\) e, portanto, temos evidência de que a escolha das pessoas pelas marcas de chocolate não é homogênea. Observando a tabela, vemos que o desequilíbrio ocorreu pela maior preferência pelo chocolate B.


Os graus de liberdade em uma tabela com \(L\) linhas e \(C\) colunas são dados por \[df = (L-1) \cdot (C-1)\] Portanto, uma tabela \(2\text{ x }2\), tem somente um grau de liberdade, e uma tabela \(4\text{ x }3\) tem seis graus de liberdade.

exemplo 2: genética

As probabilidades relacionadas com \(H_0\) não precisam ser sempre as mesmas. Por exemplo, em genética, é comum precisarmos avaliar como é a herança de um genótipo a partir dos fenótipos observados.

Duas linhagens puras de plantas foram cruzadas: uma com pétalas amarelas e outra com pétalas vermelhas. A primeira geração, \(F_1\), tem pétalas cor de laranja. Então \(F_1\) é cruzada entre si, gerando \(F_2\) com 320 plantas, das quais 182 têm pétalas cor de laranja, 61 amarelas e 77 vermelhas.

Há duas possibilidades para explicar este resultado, e queremos saber qual é o caso para o fenótipo das pétalas das flores desta espécie:

  1. dominância incompleta ou codominância, com os seguintes genótipos e frequências esperadas:
    • Y/R: laranja (\(1 \over 2\))
    • Y/Y: amarela (\(1 \over 4\))
    • R/R: vermelha (\(1 \over 4\))

Neste caso, os números esperados para \(F_2\) são, respectivamente, 160, 80 e 80, seguindo as proporções 2:1:1.

  1. epistasis recessiva, com os seguintes genótipos e frequências esperadas:
    • Y/-; R/-: laranja (\(9 \over 16\))
    • y/y; R/-: amarela (\(3 \over 16\))
    • Y/-; r/r: vermelha (\(3 \over 16\))
    • y/y; r/r: vermelha (\(1 \over 16\))

Neste caso, os números esperados para \(F_2\) são, respectivamente, 180 , 60 e 80 (60+20), seguindo as proporções 9:3:(3:1).

A qual das duas hipóteses a herança da cor das flores desta espécie adere?

Implementamos uma solução em eiras.quiaderencia.R e TesteQuiQuadradoAderencia_Flores.R.

Observe o código com nova função declarada:

# eiras.quiaderencia.R
# Teste Qui-quadrado de aderencia

eiras.quiaderencia <- function (df_dados, B=1e6)
{
  valores <- as.vector(unlist(df_dados[,1]))
  probabilidades <- as.vector(unlist(df_dados[,2]))
  # teste robusto (bootstrapping)
  robusto <- chisq.test(valores, p=probabilidades, 
                        simulate.p.value = TRUE, B = B) 
  # usando o teste convencional para obter os graus de liberdade
  tradicional <- chisq.test(valores, p=probabilidades) 
  df <- tradicional$parameter
  X2 <- robusto$statistic
  p <- robusto$p.value
  n <- sum(robusto$observed)
  V <- sqrt((X2/n)/df)
  # resultados  
  cat("----------------------------------------------------\n")
  cat("Teste de aderencia a distribuicao multinomial\n")
  prmatrix(df_dados, rowlab = rep("",nrow(df_dados)))
  cat("X^2 =", X2 , "gl =",  df, "p exato =", p, "\n")
  cat("\nHeuristica de significancia: X^2/gl > 2:\n", X2/df, "\n")
  cat("\nResiduos estandardizados:\n", robusto$stdres, "\n")
  if (0 <= V & V < 0.1) {gV <- "minimo"}
  if (0.1 <= V & V < 0.3) {gV <- "pequeno"}
  if (0.3 <= V & V < 0.5) {gV <- "intermediario"}
  if (0.5 <= V & V <= 1.0) {gV <- "grande"}
  cat("V de Cramer = ", V,"\n",sep="")
  cat("(Grau ", gV, " de dessemelhanca entre as distribuicoes observada e hipotetizada\n",sep="")
  cat("----------------------------------------------------\n")
}

Criar uma função e armazená-la em um arquivo separado facilita seu uso (o nome do arquivo tem o mesmo nome da função só por conveniência; poderia ter qualquer nome). Esta função, eiras.quiaderencia(), recebe um data frame com duas colunas e o número de iterações para o bootstrapping (com \(10^6\) por default).

Adaptar o código para outros casos torna-se muito mais simples. No código anterior, seria preciso alterar todas as ocorrências de df_chocolate e os nomes das colunas Pessoas e probH0. Agora, só é necessário modificar TesteQuiQuadradoAderencia_Flores.R:

# TesteQuiQuadradoAderencia_Flores.R
library(readxl)

source ("eiras.quiaderencia.R")

# pega o arquivo e mostra seu conteúdo
df_flores <- read_excel("flores.xlsx")
cat("\nDados:\n")
prmatrix(df_flores, rowlab = rep("",nrow(df_flores)))

# chama os testes estatísticos
cat("\nHipotese de codominancia:\n")
eiras.quiaderencia(df_flores[,c(2,3)], B=1e6)

cat("\nHipotese de epistasis recessiva:\n")
eiras.quiaderencia(df_flores[,c(2,4)], B=1e6)

Para funcionar, carrega-se a função na memória (com source()), lemos o arquivo de interesse flores.xlsx e chamamos eiras.quiaderencia() duas vezes, indicando as colunas a serem testadas:

  • a coluna 2 da planilha (que tem o número de plantas com cada característica) contra a coluna 3, que testa a hipótese de codominância,
  • a coluna 2 da planilha contra a coluna 4, que testa a hipótese de epistasis.

Este código produz a seguinte saída:


Dados:
 CorPetalas Numero Codominancia Epistasis
 "laranja"  "182"  "0.50"       "0.5625" 
 "vermelha" " 77"  "0.25"       "0.2500" 
 "amarela"  " 61"  "0.25"       "0.1875" 

Hipotese de codominancia:
----------------------------------------------------
Teste de aderencia a distribuicao multinomial
 Numero Codominancia
    182         0.50
     77         0.25
     61         0.25
X^2 = 7.65 gl = 2 p exato = 0.02224998 

Heuristica de significancia: X^2/gl > 2:
 3.825 

Residuos estandardizados:
 2.459675 -0.3872983 -2.452889 
V de Cramer = 0.1093303
(Grau pequeno de dessemelhanca entre as distribuicoes observada e hipotetizada
----------------------------------------------------

Hipotese de epistasis recessiva:
----------------------------------------------------
Teste de aderencia a distribuicao multinomial
 Numero Epistasis
    182    0.5625
     77    0.2500
     61    0.1875
X^2 = 0.1513889 gl = 2 p exato = 0.9324391 

Heuristica de significancia: X^2/gl > 2:
 0.07569444 

Residuos estandardizados:
 0.2253745 -0.3872983 0.143223 
V de Cramer = 0.01538002
(Grau minimo de dessemelhanca entre as distribuicoes observada e hipotetizada
----------------------------------------------------

Neste exemplo buscamos não rejeitar \(H_0\). Para o habitual \(\alpha=0.05\) estes resultados sugerem que a herança das cores das pétalas desta espécie é explicada por epistasis recessiva.

Teste Qui-quadrado de Pearson para Independência

O teste Qui-quadrado de Pearson permite que se descubra se existe um relacionamento ou efeito de interação entre duas variáveis qualitativas nominais com duas ou mais categorias.

Cada observação da unidade observacional deve ser alocada em apenas uma categoria de cada uma das duas variáveis nominais.

SIEGEL, S. & CASTELLAN Jr., N.J. (1988) Nonparametric statistics for the behavioral sciences. 2ª ed. NY: McGraw-Hill, p. 111

Exemplo: capacete e trauma

Considere um estudo que investiga a efetividade dos capacetes de segurança de bicicleta na prevenção de traumas cranianos. Os dados consistem de uma amostra aleatória de 793 indivíduos envolvidos em acidentes ciclísticos durante um certo período. Deseja-se testar, com \(\alpha=0.05\), se o uso do capacete tem funcionado como um fator de proteção.

Formulamos as hipóteses nula e alternativa:

\(H_0: \text{não existe associação entre uso do capacete e trauma craniano}\)

\(H_1: \text{existe associação entre uso do capacete e trauma craniano}\)

\(\alpha = 5\%\)

Os dados são:

tabela <- as.table(matrix(c(17, 138, 130, 508), nrow = 2, byrow = TRUE))
colnames(tabela) <- c("Trauma +","Trauma -")
rownames(tabela) <- c("Capacete +","Capacete -")
print (tabela)
           Trauma + Trauma -
Capacete +       17      138
Capacete -      130      508

Sob a hipótese nula (\(H_0\)), i.e., se não existe associação entre o uso de capacete e o trauma craniano nos acidentes, os valores esperados são:

chi2 <- chisq.test(tabela)
print(round(chi2$expected,3))
           Trauma + Trauma -
Capacete +   28.733  126.267
Capacete -  118.267  519.733

É possível transformar em proporções, obtendo:

tabela <- as.table(matrix(c(17, 138, 130, 508), nrow = 2, byrow = TRUE))
colnames(tabela) <- c("Trauma +","Trauma -")
rownames(tabela) <- c("Capacete +","Capacete -")
tabela_prop <- tabela/sum(tabela)
print (round(tabela_prop,3))
           Trauma + Trauma -
Capacete +    0.021    0.174
Capacete -    0.164    0.641

e, sob \(H_0\), os valores esperados são:

soma_col1 <- sum(tabela_prop[,1])
soma_col2 <- sum(tabela_prop[,2])
soma_lin1 <- sum(tabela_prop[1,])
soma_lin2 <- sum(tabela_prop[2,])
tabela_propesp <- tabela_prop
tabela_propesp[1,1] <- soma_lin1*soma_col1
tabela_propesp[1,2] <- soma_lin1*soma_col2
tabela_propesp[2,1] <- soma_lin2*soma_col1
tabela_propesp[2,2] <- soma_lin2*soma_col2
print(round(tabela_propesp,3))
           Trauma + Trauma -
Capacete +    0.036    0.159
Capacete -    0.149    0.655

Note que o teste em R é feito diretamente com os valores observados (pois o tamanho da amostra importa) e seus resultados foram armazenados, acima, na variável chi2. O tipo desta variável é dado por:

str(chi2)
List of 9
 $ statistic: Named num 6.7
  ..- attr(*, "names")= chr "X-squared"
 $ parameter: Named int 1
  ..- attr(*, "names")= chr "df"
 $ p.value  : num 0.00964
 $ method   : chr "Pearson's Chi-squared test with Yates' continuity correction"
 $ data.name: chr "tabela"
 $ observed : 'table' num [1:2, 1:2] 17 130 138 508
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:2] "Capacete +" "Capacete -"
  .. ..$ : chr [1:2] "Trauma +" "Trauma -"
 $ expected : num [1:2, 1:2] 28.7 118.3 126.3 519.7
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:2] "Capacete +" "Capacete -"
  .. ..$ : chr [1:2] "Trauma +" "Trauma -"
 $ residuals: 'table' num [1:2, 1:2] -2.189 1.079 1.044 -0.515
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:2] "Capacete +" "Capacete -"
  .. ..$ : chr [1:2] "Trauma +" "Trauma -"
 $ stdres   : 'table' num [1:2, 1:2] -2.7 2.7 2.7 -2.7
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:2] "Capacete +" "Capacete -"
  .. ..$ : chr [1:2] "Trauma +" "Trauma -"
 - attr(*, "class")= chr "htest"
sapply(chi2, typeof)
  statistic   parameter     p.value      method   data.name    observed 
   "double"   "integer"    "double" "character" "character"    "double" 
   expected   residuals      stdres 
   "double"    "double"    "double" 

É uma estrutura complexa, da qual muitos valores podem ser extraídos quando quisermos utilizá-los dentro de um RScript. A parte em chi2$expected, que já utilizamos, contém a tabela com os valores esperados sob \(H_0\).

O resultado principal do teste é exibido pela variável como um todo:

print (chi2)

    Pearson's Chi-squared test with Yates' continuity correction

data:  tabela
X-squared = 6.7001, df = 1, p-value = 0.009641

A estatística de teste tem o valor \(X^2=6.7001\), com 1 grau de liberdade (\(df=1\)), correspondendo ao valor-p de \(0.009641\) (que será discutido adiante). Com estes valores podemos decidir pela rejeição ou não-rejeição de \(H_0\).

No exemplo de uma tabela 2x2 há somente 1 grau de liberdade. O gráfico da distribuição de qui-quadrado (para melhor visibilidade, computamos \(0 \le \chi^2 \le 8\), pois o valor de \(\chi^2\) vai de zero a \(\infty\), e truncamos o eixo \(y\) em 1.2) é:

chi2.valor <- seq(from=0, to=8, by=0.01)
graus.liberdade <- 1
p <- dchisq(chi2.valor, df=graus.liberdade)
plot(chi2.valor, p, type="l", ylim=c(0, 1.2),
     xlab="Qui-quadrado", ylab="Densidade"
)


Conceitualmente, o teste Qui-quadrado é um teste de proporções, no qual a estatística de teste (i.e., o valor de \(X^2\)) é obtida pela diferença quadrática entre as frequências absolutas observadas e esperadas, dada por \[X^2 = \sum{ {(o - e)^2} \over{e} }\]. ou \[X^2 = {{(ad-bd)^2(a+b+c+d)} \over {(a+b)(c+d)(b+d)(a+c)}}\]

Em R, experimente calcular a estatística na “força bruta”:

qui2 <- 0
for(i in 1:4)
{
  qui2 <- qui2 + ((chi2$observed[i]-chi2$expected[i])^2)/chi2$expected[i]
}
cat("O valor do qui-quadrado é ",qui2,"\n",sep="")
O valor do qui-quadrado é 7.309882

O valor de \(X^2\) aumenta com a discrepância crescente entre os valores observados e esperados.

Note, porém, que o valor calculado pelo código acima é diferente do que aparece em:

print (chi2)

    Pearson's Chi-squared test with Yates' continuity correction

data:  tabela
X-squared = 6.7001, df = 1, p-value = 0.009641

Isto acontece por causa da correção para continuidade de Yates, que será comentada adiante. Na “força bruta”, a correção pode ser implementada com:

qui2y <- 0
for(i in 1:4)
{
  qui2y <- qui2y + ((abs(chi2$observed[i]-chi2$expected[i])-0.5)^2)/chi2$expected[i]
}
cat("O valor do qui-quadrado com a correção de Yates é ",qui2y,"\n",sep="")
O valor do qui-quadrado com a correção de Yates é 6.700121
Obviamente, nenhum destes cálculos manuais é necessário. A função chisq.test() já os executa e armazena no objeto que denominamos, neste exemplo, chi2.

Nosso interesse é saber a partir de qual o valor de \(X^2\) a discrepância entre os valores observados e esperados é menos provável que o valor de \(\alpha\) escolhido. Isto é encontrar o valor que separa a área de \(5\%\) sob a curva de distribuição em sua cauda direita. Este valor é chamado de qui-quadrado crítico, \(\chi^2_c\).


Note que o raciocínio é bicaudal, buscando associação positiva ou negativa, mas a operacionalização do teste só utiliza a cauda superior. Isto acontece porque \(chi^2=0\) indica igualdade numérica (e, portanto, também estatística) entre os valores esperados e observados. Consequentemente, procuramos \(\alpha=0.05\) somente na cauda direita.

Em R, podemos encontrar \(\chi^2_c\) pelo complemento (\(1-\alpha=0.95\)) com

chi2.critico <- qchisq(p=0.95, df=1)
cat("qui-quadrado crítico (df=1, alfa=5%): ", chi2.critico, "\n", sep="")
qui-quadrado crítico (df=1, alfa=5%): 3.841459

pois este é o valor de \(\chi^2\) que deixa \(95\%\) da área sob a curva à sua esquerda.

Graficamente, podemos hachurar em azul a área correspondente a \(\alpha\) e traçar uma linha pontilhada vertical em \(\chi^2_c\) com:

source("friendlycolor.R")
# vetor com valores de 0 a 8
chi2.valor <- seq(from=0, to=8, by=0.01)
graus.liberdade <- 1
# vetor com probabilidades correspondentes
p <- dchisq(chi2.valor, df=graus.liberdade)
# exibe o gráfico
plot(chi2.valor, p, type="l", ylim=c(0, 1.2),
     xlab="Qui-quadrado", ylab="Densidade")
# vetor com os valores acima de qui critico 
chi2.critico <- qchisq(p=0.95, df=1)
chi2maioralfa <- chi2.valor[chi2.valor>chi2.critico] 
# probabilidades correspondentes
pmaioralfa <- dchisq(chi2maioralfa, df=graus.liberdade)
# hachura sob a curva
chi2maioralfa <- c(min(chi2maioralfa), chi2maioralfa, max(chi2maioralfa))
pmaioralfa <- c(0, pmaioralfa, 0)
polygon(chi2maioralfa, pmaioralfa, col=friendlycolor(10), border = NA)
# reforca a linha
lines(chi2maioralfa[2:(length(chi2maioralfa)-1)], 
      pmaioralfa[2:(length(pmaioralfa)-1)], 
      col=friendlycolor(10), lwd=2)
# linha pontilhada, valor critico
abline(v=chi2.critico, lty=2)

Esta área hachurada é a região de rejeição de \(H_0\).

- decisão pela estatística de teste

A estatística \(\chi^2\) é maior quanto maior for a discrepância entre os valores observados e esperados. Para um grau de liberdade e \(\alpha=0.05\), computamos \(\chi^2_c=\) 3.8414588. No exemplo de capacete e trauma, acima, rejeitamos \(H_0\) porque calculamos \(X^2=\) 6.7001206, valor maior que o valor crítico.

- decisão pelo valor-p

À direita do valor calculado \(X^2=\) 6.7001206 define-se uma área, que corresponde ao valor-p. Este valor também foi calculado quando aplicamos o teste e guardamos seu resultado em chi2, acima; o valor-p encontra-se acessível em chi2$p.value=0.0096406. Como \(p < \alpha\), concluimos pela rejeição de \(H_0\).

Podemos hachurar em laranja a área correspondente ao valor-\(p\) e traçar uma linha pontilhada vertical na altura do \(X^2\) calculado, adicione mais linhas ao código anterior, mas precisamos da variável chi2$statistic que retorna de chisq.tes(). O gráfico, agora, tem muitos elementos e, portanto, também devemos adicionar uma legenda.

O código completo é:

tabela <- as.table(matrix(c(17, 138, 130, 508), nrow = 2, byrow = TRUE))
colnames(tabela) <- c("Trauma +","Trauma -")
rownames(tabela) <- c("Capacete +","Capacete -")
chi2 <- chisq.test(tabela)
source("friendlycolor.R")
# vetor com valores de 0 a 8
chi2.valor <- seq(from=0, to=8, by=0.01)
graus.liberdade <- 1
# vetor com probabilidades correspondentes
p <- dchisq(chi2.valor, df=graus.liberdade)
# exibe o gráfico
plot(chi2.valor, p, type="l", ylim=c(0, 1.2),
     xlab="Qui-quadrado", ylab="Densidade")
# vetor com os valores acima de qui critico 
chi2.critico <- qchisq(p=0.95, df=1)
chi2maioralfa <- chi2.valor[chi2.valor>chi2.critico] 
# probabilidades correspondentes
pmaioralfa <- dchisq(chi2maioralfa, df=graus.liberdade)
# hachura sob a curva
chi2maioralfa <- c(min(chi2maioralfa), chi2maioralfa, max(chi2maioralfa))
pmaioralfa <- c(0, pmaioralfa, 0)
polygon(chi2maioralfa, pmaioralfa, col=friendlycolor(10), border = NA)
# reforca a linha
lines(chi2maioralfa[2:(length(chi2maioralfa)-1)], 
      pmaioralfa[2:(length(pmaioralfa)-1)], 
      col=friendlycolor(10), lwd=2)
# linha pontilhada, valor critico
abline(v=chi2.critico, lty=2)

# ########## Area a direita de p ############
# vetor com os valores acima de qui calculado 
chi2maiorp <- chi2.valor[chi2.valor>chi2$statistic] 
# probabilidades correspondentes
pmaiorp <- dchisq(chi2maiorp, df=graus.liberdade)
# hachura sob a curva
chi2maiorp <- c(min(chi2maiorp), chi2maiorp, max(chi2maiorp))
pmaiorp <- c(0, pmaiorp, 0)
polygon(chi2maiorp, pmaiorp, col=friendlycolor(20), border = NA)
# reforca a linha
lines(chi2maiorp[2:(length(chi2maiorp)-1)], 
      pmaiorp[2:(length(pmaiorp)-1)], 
      col=friendlycolor(20), lwd=2)
# linha pontilhada, valor critico
abline(v=chi2$statistic, lty=3)

# legenda
legend("topright",
       c("Qui^2 crítico","Qui^2 obs.", "alfa", "p"), 
       col=c("black", "black", friendlycolor(10), friendlycolor(20)),
       lty=c(2, 3, 1, 1), 
       lwd=c(1, 1, 5, 5),
       box.lwd=0, bg="white",
       cex=0.8)  


As duas formas de decisão são igualmente válidas e equivalentes.

Exibir o valor-p é a forma mais moderna e recomendada atualmente.

Teste Qui-quadrado robusto

Repare nas saídas da função chisq.test(), cujo cabeçalho indica:

Pearson’s Chi-squared test with Yates’ continuity correction

Esta é uma versão tradicional do teste, para a qual a correção de Yates é aplicável por default às tabelas 2x2 (mais detalhes serão vistos adiante). Podemos evitá-la com:

chi2ny <- chisq.test(tabela, correct=FALSE)
print(chi2ny)

    Pearson's Chi-squared test

data:  tabela
X-squared = 7.3099, df = 1, p-value = 0.006858

Além disto, há várias condições para a aplicação do teste Qui-quadrado (comentadas adiante) que podem, muitas vezes, restringir seu uso.

Uma alternativa é computar o teste em sua versão robusta. A função chisq.test() consegue simular a distribuição por bootstrapping, alterando-se o comando para:

chi2.robusto <- chisq.test(tabela, correct=FALSE,
                   simulate.p.value = TRUE, B=1e6)
print(chi2.robusto)

    Pearson's Chi-squared test with simulated p-value (based on 1e+06
    replicates)

data:  tabela
X-squared = 7.3099, df = NA, p-value = 0.007638

Esta versão é simulada com \(1e6 = 1 \cdot 10^6 = 1000000\) iterações.


Como a amostra é grande, os valores de \(X^2\) e \(p\) não mudam muito entre as diferentes versões. Há um pequeno bug na implementação da versão robusta, que não calcula os graus de liberdade (sabemos que \(df=1\)).

A decisão, neste caso, não mudou, tanto pelo critério do \(\chi^2_c\) quanto pelo \(p\): rejeita-se \(H_0\).

Riscos relativos


Odds e probabilidade

Odds (traduzido habitualmente por chance) e probabilidade são formas equivalentes para expressar possibilidades e relacionadas por:

\(\textit{Odds} = {{\text{Probabilidade}} \over{1-\text{Probabilidade}}}\)

e

\(\text{Probabilidade} = {{\textit{Odds}} \over{1+\textit{Odds}}}\)

Por exemplo, uma probabilidade de \(80\%\) corresponde a

\(\textit{Odds} = {{0.8} \over {1-0.8}} = {0.8 \over 0.2} = 4\)

indicando que dizer “\(80\%\) de probabilidade” é equivalente a dizer “4 vezes mais chance de ocorrer do que não ocorrer” um determinado evento.

Resersamente, Odds de 2 é:

\(\text{Probabilidade} = {2 \over {1+2}} = {2 \over 3} \approx 66.67\%\)

e \(\textit{Odds}=1\) resulta em:

\(\text{Probabilidade} = {1 \over {1+1}} = {1 \over 2} = 50\%\)

Associa-se o máximo de incerteza com probabilidade de \(50\%\); por este motivo, \(Odds=1\) será um valor necessário para decisões estatísticas, adiante.

Há dois tipos principais de risco relativo: razão de riscos (RR, risk ratio) e razão de chances (OR, odds ratio).

Vamos considerar a seguinte tabela:

                  Tem efeito Não tem efeito      
Tem exposição     a          b              a+b  
Não tem exposição c          d              c+d  
                  a+c        b+d            total

Esta tabela de contingência genericamente relaciona exposição (e.g., hábito de fumar) com o efeito (e.g., câncer de pulmão). Traz as contagens dos:

  • expostos que apresentaram o efeito (\(a\)),
  • expostos que não apresentaram o efeito (\(b\)),
  • não expostos que apresentaram o efeito (\(c\)) e
  • não expostos que não apresentaram o efeito (\(d\)),

razão de riscos (RR, risk ratio)

Riscos são probabilidades. A razão de riscos é dada por:

\(RR = \LARGE{{{{a} \over {a+b}} \over {{{c} \over {c+d}}}}}\)

Repare que \({a} \over {a+b}\) é a probabilidade de um indivíduo exposto ter o efeito (e.g., fumante ter câncer); \({c} \over {c+d}\) é a probabilidade de um indivíduo não exposto ter o efeito (e.g., não-fumante ter câncer).

RR, portanto, é um odds: quantas vezes mais é esperado o efeito (ocorrência de câncer) de quem é exposto (fumante) em comparação com quem não é exposto (não fumante);

razão de chances (OR, odds ratio)

Odds ratio como diz o nome, é uma razão entre dois odds, dada por:

\(OR = \LARGE{{ { { {a} \over {a+c} } \over { {c} \over {a+c} } } \over { { {b} \over {b+d} } \over { {d} \over {b+d} } } }}\)

Considere, agora, que:

  • \(\LARGE{ { {a} \over {a+c} } \over { {c} \over {a+c} } }\), entre os doentes, quantas vezes é a maior a chance (odds) de ter efeito (câncer) para quem é exposto (fuma).
  • \(\LARGE{ { {b} \over {b+d} } \over { {d} \over {b+d} } }\), entre os não doentes, quantas vezes é a maior a chance (odds) de ter efeito (câncer) para quem é exposto (fuma).

Interessante, também, é que (com alguma álgebra) \(OR\) reduz-se a:

\(OR = \LARGE{ {{a \cdot d} \over {{b \cdot c}} } }\)

fácil de lembrar e assim conhecido como “produto cruzado” (com o defeito de ocultar o motivo pelo qual é um odds ratio). Por outro lado, como produto cruzado, percebe-se que é uma medida de quanto pesa a:

  • diagonal principal, do que é concordante: \(\text{efeito entre expostos}~~\text{E}~~\text{não efeito entre não expostos}\)

em relação à

  • diagonal secundária, do que é discordante: \(\text{efeito entre não expostos}~~\text{E}~~\text{não efeito entre expostos}\).

Os valores de \(RR\) e \(OR\) são próximos quando menos que \(5\%\) dos não expostos têm efeito.

Em R, o cálculo de \(OR\) aparece no Teste Exato de Fisher. No exemplo de capacete e trauma, temos:

tabela <- as.table(matrix(c(17, 138, 130, 508), nrow = 2, byrow = TRUE))
colnames(tabela) <- c("Trauma +","Trauma -")
rownames(tabela) <- c("Capacete +","Capacete -")
print(tabela)
ft <- fisher.test(tabela) # Teste de OR robusto
print(ft)
           Trauma + Trauma -
Capacete +       17      138
Capacete -      130      508

    Fisher's Exact Test for Count Data

data:  tabela
p-value = 0.005688
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 0.2630993 0.8357342
sample estimates:
odds ratio 
 0.4817645 

- decisão pelo intervalo de confiança

A estimativa pontual de odds ratio está em ft$estimate=0.4817645, mas para a decisão é obrigatório verificar seu intervalo de confiança 95% (\(IC95\)), com limites inferior e superior respectivamente guardados em ft$conf.int[1]=0.2630993 e ft$conf.int[2]=0.8357342.

Como a ausência de efeito é dada por \(OR=1\) (correspondendo a \(H_0\)) e este valor unitário está fora do intervalo, rejeita-se \(H_0\). Além disto, o \(IC95\) é menor que o esperado por \(H_0\) e, portanto, o uso de capacete está associado com redução dos traumas cranianos.


Caso a tabela fosse construída com as colunas em posições trocadas, obter-se-ia o seguinte:

           Trauma - Trauma +
Capacete +      138       17
Capacete -      508      130

    Fisher's Exact Test for Count Data

data:  tabela
p-value = 0.005688
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 1.196553 3.800847
sample estimates:
odds ratio 
  2.075703 

O odds ratio é agora maior que o valor unitário, o qual está fora do intervalo, rejeitando-se a hipótese nula e a associação é positiva: a conclusão é a mesma, o uso de capacete está associado com o aumento de “Trauma -”, portanto é protetor.

Aliás, observe que

print(1/ft$estimate)
odds ratio 
 0.4817645 
print(1/ft$conf.int)
[1] 0.8357342 0.2630993
attr(,"conf.level")
[1] 0.95
são os mesmos valores obtidos anteriormente.

\(RR\) e \(OR\) em epitools

Há outros pacotes em R que calculam risk ratio e odds ratio, como por exemplo:

library(epitools)

tabela <- as.table(matrix(c(17, 138, 130, 508), nrow = 2, byrow = TRUE))
colnames(tabela) <- c("Trauma +","Trauma -")
rownames(tabela) <- c("Capacete +","Capacete -")

epitools::oddsratio(tabela, method="fisher")

que produz os mesmos valores para a estimativa pontual e intervalo de confiança que o teste exato de Fisher:

$data
           Trauma + Trauma - Total
Capacete +       17      138   155
Capacete -      130      508   638
Total           147      646   793

$measure
                        NA
odds ratio with 95% C.I.  estimate     lower     upper
              Capacete + 1.0000000        NA        NA
              Capacete - 0.4817645 0.2630993 0.8357342

$p.value
            NA
two-sided     midp.exact fisher.exact  chi.square
  Capacete +          NA           NA          NA
  Capacete - 0.005088331  0.005688151 0.006857643

$correction
[1] FALSE

attr(,"method")
[1] "Conditional MLE & exact CI from 'fisher.test'"

Há mais três métodos disponíveis, de acordo com a documentação da função.

A library epitools também permite o cálculo de \(RR\) com intervalo de confiança. Ainda que não seja adequada para o presente exemplo, verifiquem sua sintaxe. Há três métodos disponíveis, como por exemplo:

library(epitools)

tabela <- as.table(matrix(c(17, 138, 130, 508), nrow = 2, byrow = TRUE))
colnames(tabela) <- c("Trauma +","Trauma -")
rownames(tabela) <- c("Capacete +","Capacete -")

epitools::riskratio(tabela, method="wald")
$data
           Trauma + Trauma - Total
Capacete +       17      138   155
Capacete -      130      508   638
Total           147      646   793

$measure
                        NA
risk ratio with 95% C.I.  estimate     lower     upper
              Capacete + 1.0000000        NA        NA
              Capacete - 0.8943256 0.8357184 0.9570428

$p.value
            NA
two-sided     midp.exact fisher.exact  chi.square
  Capacete +          NA           NA          NA
  Capacete - 0.005088331  0.005688151 0.006857643

$correction
[1] FALSE

attr(,"method")
[1] "Unconditional MLE & normal approximation (Wald) CI"

e também uma variante com bootstrapping:

library(epitools)

tabela <- as.table(matrix(c(17, 138, 130, 508), nrow = 2, byrow = TRUE))
colnames(tabela) <- c("Trauma +","Trauma -")
rownames(tabela) <- c("Capacete +","Capacete -")

epitools::riskratio(tabela, replicates=1e6)
$data
           Trauma + Trauma - Total
Capacete +       17      138   155
Capacete -      130      508   638
Total           147      646   793

$measure
                        NA
risk ratio with 95% C.I.  estimate     lower     upper
              Capacete + 1.0000000        NA        NA
              Capacete - 0.8943256 0.8357184 0.9570428

$p.value
            NA
two-sided     midp.exact fisher.exact  chi.square
  Capacete +          NA           NA          NA
  Capacete - 0.005088331  0.005688151 0.006857643

$correction
[1] FALSE

attr(,"method")
[1] "Unconditional MLE & normal approximation (Wald) CI"

Lembre-se de nunca decidir a respeito de risk ratio ou odds ratio sem considerar o intervalo de confiança. Por isso, evite fazer o cálculo manualmente: use R!

Dois códigos úteis, com vários dos elementos apresentados neste texto

4
Simulação e Gráficos

1. tabelacontingencia2x2_capacetetrauma.R

Em tabelacontingencia2x2_capacetetrauma.R, que você pode adaptar para as suas necessidades, processamos o mesmo exemplo (procure a função sink() que desvia o resultado da tela para um arquivo texto). Além da estatística \(X^2\), computamos outros valores interessantes. Observe sua saída (tabelacontingencia2x2_capacetetrauma.txt):

source("tabelacontingencia2x2_capacetetrauma.R")

Resultados armazenados em tabelacontingencia2x2_capacetetrauma.txt

Em Qui2.R também processamos o mesmo exemplo. Esta versão aplica bootstrapping de um jeito mais manual e produz gráficos associados aos conceitos que discutimos.

Execute com:

source("Qui2.R")
fornecendo:
**** TABELA DE CONTINGÊNCIA ****

Considere:
                                               
              Coluna+      Coluna-             
 Linha +            a            b Total lin. +
 Linha -            c            d Total lin. -
         Total col. + Total col. -  Total geral


Informe:
Nome nas linhas (ate 10 letras): Capacete
Nome nas colunas (ate 10 letras): Trauma
Definido:
                                                                     
                       Trauma +            Trauma -                  
 Capacete +                   a                   b Total de Trauma +
 Capacete -                   c                   d Total de Trauma -
            Total de Capacete + Total de Capacete -       Total geral

Forneca os valores de a, b, c, d:
a: 17
b: 138
c: 130
d: 508

Defina alfa = prob. erro do tipo I).
(número entre 0 e 1, default = 0.05).
alfa: 0.05

Quantas iteracoes para simular?
(entre um numero inteiro; default n=1e4)
iteracoes: 1e4

Exibir tabelas? (exibir lentifica a simulacao)
0=nao, 1=sim; default eh 0: 0

Observe a simulação em andamento na área de Plots. Quando terminar, procure os gráficos finalizados. Os principais são:

  • as distribuições de \(\chi^2\) para a hipótese nula e para o valor observado de \(H_1\):

O valor de \(\beta\) e o poder associados, aqui, são os valores a posteriori, portanto inúteis para qualquer decisão


Flashback:

Recorde a tomada de decisão com a distribuição binomial: A estatística e o formato das distribuições mudam (em função, também, do tipo de variável), mas as áreas de \(\alpha\), \(\beta\) e o raciocínio são os mesmos.

  • a distribuição, com intervalo de confiança, do odds ratio simulado (compare com o obtido acima), indicando a rejeição de \(H_0\) (o valor 1, de referência, está fora do intervalo) e que seu efeito está entre mínimo a intermediário:

Odds ratio é uma medida da intensidade de associação entre duas variáveis nominais dicotômicas de exposição e desfecho.

De acordo com SHARPE, D. (2015) Your chi-square test is statistically significant: now what? Practical Assessment, Research & Evaluation, 20(8). É, também, uma medida de tamanho de efeito (não depende do tamanho da amostra, é adimensional e adirecional) e sua intensidade pode ser classificada em:

  • OR para categoria de proteção
    • Grande: \(OR < 0,1\)
    • Intermediária: \(0,1 \le OR < 0,3\)
    • Pequena: \(0,3 \le OR < 0,7\)
    • Mínima: \(0,7 \le OR < 1\)
  • OR para categoria de risco
    • Mínima: \(1 < OR \le 1,5\)
    • Pequena: \(1,5 < OR \le 3,5\)
    • Intermediária: \(3,5 < OR \le 9\)
    • Grande: \(OR > 9\)

  • a distribuição do tamanho de efeito (V de Cramer), mostrando que o efeito, apesar de significante, está entre mínimo e pequeno:

O programa já indica o tamanho de efeito de acordo com O V de Cramer é uma medida do grau de correlação absoluta entre duas variáveis nominais:

  • não depende do total de observações independentes
  • não depende do número de linhas e colunas da TC
  • é adimensional e varia entre 0 (independência) e 1 (dependência)

Tamanho de efeito:

  • \(0 \le V < 0.1\) : mínimo
  • \(0.1 \le V < 0.3\) : pequeno
  • \(0.3 \le V < 0.5\) : intermediário
  • \(V \ge 0.5\) : grande

  • novamente a distribuição de \(\chi^2\), mas aqui exibe somente \(H_0\) e o valor-p, indicando a rejeição de \(H_0\):

2. tabelas LxC

O teste Qui-quadrado opera em tabelas maiores, com \(L\) linhas e \(C\) colunas.

Por exemplo, três quimioterápicos foram comparados quanto ao efeito colateral (náusea). Considere \(\alpha=0.05\). Há diferença entre as drogas?

Os dados obtidos foram:
            Nausea   Sem_nausea
Droga A        3         5
Droga B        7         2
Droga C        6         3

Em Qui2_LxC.R implementamos este exemplo (versão robusta):

source("Qui2_LxC.R")
  ComNausea SemNausea
A         3         5
B         7         2
C         6         3

    Pearson's Chi-squared test with simulated p-value (based on 1e+06
    replicates)

data:  tcrxc
X-squared = 3.0559, df = NA, p-value = 0.2744

V de Cramer = 0.3428334

Observe que não rejeita-se \(H_0\). Os dois critérios são sempre equivalentes:

  • o \(X^2\) calculado é menor que \(\chi^2_c\)
  • \(p > \alpha\)

Abra o código para ver como foi implementado este teste.

Comentários sobre os métodos tradicionais

5
Condições para o teste Qui-quadrado
(métodos não robustos)

Antes das implementações computacionais e a facilidade do uso do R, a partir de uma tabela 2x2 era necessário construir a tabela dos valores esperados.

  • Valores observados:
           Trauma + Trauma -
    Capacete +       17      138     155
    Capacete -      130      508     638
                147      646     793
    
  • Valores esperados:
     
           Trauma + Trauma -
    Capacete +     28.7    126.3
    Capacete -    118.3    519.7
    

    Cada valor esperado é obtido pelo produto das marginais, dividido pelo total geral:

    • \((147 \cdot 155) / 793 \approx 28.7\)
    • \((147 \cdot 638) / 793 \approx 118.3\)
    • \((646 \cdot 155) / 793 \approx 126.3\)
    • \((646 \cdot 638) / 793 \approx 519.7\)

Avaliando os dados observados e os valores esperados, era necessário verificar se o teste Qui-quadrado podia ser aplicado. Eram as seguintes condições:

SIEGEL, S. & CASTELLAN Jr., N.J. (1988) Nonparametric statistics for the behavioral sciences. 2ª ed. NY: McGraw-Hill, p. 123:

  • O teste exato de Fisher para testar independência é adequado apenas para tabela de contingência 2×2 se N < 20
  • Se N entre 20 e 40, o teste Qui-quadrado de Pearson pode ser usado se todas as frequências esperadas são maiores que 5.
  • Se N > 40, o teste Qui-quadrado de Pearson com correção de Yates pode ser usado.

COCHAN, W.G. (1954) Some Methods for Strengthening the Common \(\chi^2\) Tests. Biometrics: 10(4): 417-451:

  • Um valor esperado mínimo de 1 em alguma célula é permitido, desde que não mais que 20% das células tenham valor abaixo de 5

Então, o valor de \(X^2\) era dado por

\[X^2 = \sum{ ({\text{observado}-\text{esperado}})^2 \over {\text{observado}} }\]

neste exemplo calculado por:

\(X^2 = { {{(17-28.7)^2} \over {28.7}} + {{(138-126.3)^2} \over {126.3}} + {{(130-118.3)^2} \over {118.3}} + {{(508-519.7)^2} \over {519.7}} } \approx 7.31\)

ou, em tabelas 2x2, era recomendado usar a correção de continuidade de Yates, dada por:

\[X^2 = \sum{ (|{\text{observado}-\text{esperado}}|-0.5)^2 \over {\text{observado}} }\]

Neste exemplo:

\(X^2 = { {{(|17-28.7|-0.5)^2} \over {28.7}} + {{(|138-126.3|-0.5)^2} \over {126.3}} + {{(|130-118.3|-0.5)^2} \over {118.3}} + {{(|508-519.7|-0.5)^2} \over {519.7}} } \approx 6.7\)

Então, era necessário consultar uma tabela com os valores críticos de \(\chi^2_c\) previamente calculados para a tomada de decisão (sem obter, portanto, o valor-\(p\)). A tabela permitia escolher o valor de \(\alpha\) e, sabendo-se os graus de liberdade (neste caso \(\nu = 1\)), localizar \(\chi^2_c\) (neste exmeplo igual a \(3.841\)):

Caso \(\chi^2\) não atendesse as exigências para ser aplicado, a alternativa (se fosse tabela 2x2) era o Teste Exato de Fisher. Este envolve calcular o valor-\(p\) por fatoriais de tabelas progressivamente mais extremas. Por exemplo: Neste caso o valor exato de \(p\) era calculado e comparado diretamente com \(\alpha\) para a tomada de decisão.

O que diz a literatura mais recente sobre estes métodos tradicionais e não robustos?

Métodos avançados

Análise post hoc

Esta análise permite, quando se rejeita \(H_0\), localizar quais células da tabela de contingência mais contribuíram para esta rejeição.

exemplo: etilismo e tabagismo

Para descobrir se existe uma relação entre tabagismo e etilismo a partir de um estudo realizado com 100 estudantes universitários, encontrou-se:

Tabagista Não tabagista
Etilista 50 15
Não etilista 20 25

A manipulação de dados em R é poderosa, mas nem sempre fácil. Para a entrada de dados na forma de uma matriz 2x2 o seguinte código funciona:

Tabela <- ("
 TC2x2       Tabagista NaoTabagista
 Etilista    50        15
 NaoEtilista 20        25  
")
TC <- as.matrix(read.table(textConnection(Tabela), 
                               header=TRUE, row.names=1))
print(TC)

Obtendo-se:

            Tabagista NaoTabagista
Etilista           50           15
NaoEtilista        20           25

O problema está na entrada de dados. Somente espaços em branco podem ser usados para a montagem da tabela. O alinhamento é visual, mas não necessário para o R. Funciona igualmente

Tabela <- ("
Tabagista NaoTabagista
Etilista 50 15
NaoEtilista 20  25  
")
TC <- as.matrix(read.table(textConnection(Tabela), 
                                 header=TRUE, row.names=1))
print(TC)
            Tabagista NaoTabagista
Etilista           50           15
NaoEtilista        20           25

Portanto, pode ser ruim a visualização de sua entrada.

Uma opção mais “limpa” é sempre colocar os dados em planilhas no formato .xls ou .xlsx e usar a função readxl::read_excel() para utilizá-los em R: Esta planilha é tabagismo_e_etilismo_2x2.xlsx, lida por

TCnew <- read_excel("tabagismo_e_etilismo_2x2.xlsx")
New names:
* `` -> ...1
print(TCnew)
# A tibble: 2 x 3
  ...1        Tabagista NaoTabagista
  <chr>           <dbl>        <dbl>
1 Etilista           50           15
2 NaoEtilista        20           25

A função read_excel() avisa que teve que suprir um nome (por causa da célula A1 que está vazia), mas esta célula será desprezada adiante.

No entanto, temos um problema a resolver. A função read_excel() produz um data frame e a função chisq.test() precisa de uma matriz para sua entrada:

is.data.frame(TCnew)
[1] TRUE
is.matrix(TCnew)
[1] FALSE

Temos, portanto, que converter TCnew para a representação em matrix. Existe uma função para isto:

TCmatrix <- data.matrix(TCnew)
Warning in data.matrix(TCnew): NAs introduced by coercion
is.matrix(TCmatrix)
[1] TRUE
print(TCmatrix)
     ...1 Tabagista NaoTabagista
[1,]   NA        50           15
[2,]   NA        20           25

Ainda não é o desejado. TCmatrix tem 2 linhas mas 3 colunas, além de ter perdido os nomes das linhas que TCnew trazia em sua primeira coluna:

nrow(TCmatrix)
[1] 2
ncol(TCmatrix)
[1] 3
rownames(TCmatrix)
NULL
print(TCnew[,1])
# A tibble: 2 x 1
  ...1       
  <chr>      
1 Etilista   
2 NaoEtilista

A solução, portanto, é transferir o conteúdo de TCnew como nomes das linhas de TCmatrix e, então, deletar a coluna 1 de TCmatrix que contém valores NA.

rownames(TCmatrix) <-as.character(unlist(TCnew[1:2,1]))
TCmatrix <- TCmatrix[,-1]
print(TCmatrix)
            Tabagista NaoTabagista
Etilista           50           15
NaoEtilista        20           25

Colocando tudo junto, para terminar com uma matriz chamada TC (como era originalmente), podemos ler a planilha em uma variável (e.g., TCtmp), reorganizar seu conteúdo em TC e remover TCtmp da memória com a função remove(). O código enxuto é:

TCtmp <- read_excel("tabagismo_e_etilismo_2x2.xlsx")
New names:
* `` -> ...1
TC <- data.matrix(TCtmp)
Warning in data.matrix(TCtmp): NAs introduced by coercion
rownames(TC) <-as.character(unlist(TCtmp[1:2,1]))
TC <- TC[,-1]
remove(TCtmp)
print(TC)
            Tabagista NaoTabagista
Etilista           50           15
NaoEtilista        20           25

Uma mensagem:

Warning in data.matrix(TCnew): NAs introduced by coercion

apareceu aqui. Warnings são avisos do R que não impedem o funcionamento do programa (avalie, caso a caso, se precisa corrigir algo; aqui não é necessário). Adiante mostramos como é removê-la.


Neste exemplo, o objetivo é executar o teste Qui-quadrado mas, no caso de rejeitar-se \(H_0\), executar uma análise post hoc. Observe o código de TabelaContingencia2x2_TabagismoEtilismo.R que pega os dados da planilha tabagismo_e_etilismo_2x2.xlsx:

# TabelaContingencia2x2_TabagismoEtilismo.R

# desabilita warnings
options(warn=-1)

# Tabela de Dancey&Reidy (2019, p. 275)
TCtmp <- read_excel("tabagismo_e_etilismo_2x2.xlsx")
TC <- data.matrix(TCtmp)
rownames(TC) <-as.character(unlist(TCtmp[1:2,1]))
TC <- TC[,-1]
remove(TCtmp)
cat("\n------------------------------------------\n")
cat("Dados:")
cat("\n------------------------------------------\n")
print(TC)

cat("\n------------------------------------------\n")
cat("Analise de significancia estatistica:")
cat("\n------------------------------------------\n")
alfa <- 0.05

cat("\nTeste qui-quadrado de Pearson exato:\n")
print(res <- chisq.test(TC,simulate.p.value=TRUE,B=1e6))
cat("X^2 critico de 95% = ",qchisq(p=1-alfa,df=1), "\n", sep="")
N <- sum(res$observed)
nL <- nrow(TC) # ou dim(TC)[1]
nC <- ncol(TC) # ou dim(TC)[2]
df <- (nL-1)*(nC-1)
X2 <- res$statistic # estatistica de teste qui-quadrado
cat("Graus de liberdade (nao fornecidos por bootstrapping): ", df, "\n", sep="")
cat("Heuristica de significancia (rej. H0 se X^2/gl > 2): ", X2/df, "\n", sep="")

cat("\n------------------------------------------\n")
cat("Analise post hoc:")
cat("\n------------------------------------------\n")
cat("\nResiduos ajustados standardizados corrigidos por momento (MCSTARs):\n")
STAR <- res$stdres # STandardized Adjusted Residual (STAR)
MCSTAR <- STAR/(sqrt((1-1/nL)*(1-1/nC))) # Moment-correct (MCSTAR)
print(MCSTAR)
alfaBonf <- alfa/df
zcrit <- abs(qnorm(alfaBonf/2))
cat("\n|MCSTAR critico| (alfaBonferroni=5%/",df,") = ",zcrit,"\n", sep="")

cat("\n------------------------------------------\n")
cat("Analise de significancia pratica:")
cat("\n------------------------------------------\n")
phi2 <- X2/N # phi ou w de Cohen
Dim <- min(nL,nC)-1 # dimensoes da TC: dim = max(phi^2)
V <- sqrt(phi2/Dim) # V ou C de Cramer
if (0 <= V & V < 0.1) {gV <- "minimo"}
if (0.1 <= V & V < 0.3) {gV <- "pequeno"}
if (0.3 <= V & V < 0.5) {gV <- "intermediario"}
if (0.5 <= V & V <= 1.0) {gV <- "grande"}
cat("\nV de Cramer =", V, "\nGrau", gV, 
    "de dependencia entre as duas variaveis nominais\n")

cat("\nTeste de razao de chances (OR) robusto:\n")
resft <- fisher.test(TC,conf.level = 1-alfa) # Teste de OR robusto
print(resft)

# reabilita warnings
options(warn=0)

que produz a seguinte saída:

New names:
* `` -> ...1

------------------------------------------
Dados:
------------------------------------------
            Tabagista NaoTabagista
Etilista           50           15
NaoEtilista        20           25

------------------------------------------
Analise de significancia estatistica:
------------------------------------------

Teste qui-quadrado de Pearson exato:

    Pearson's Chi-squared test with simulated p-value (based on 1e+06
    replicates)

data:  TC
X-squared = 12.121, df = NA, p-value = 0.000666

X^2 critico de 95% = 3.841459
Graus de liberdade (nao fornecidos por bootstrapping): 1
Heuristica de significancia (rej. H0 se X^2/gl > 2): 12.12149

------------------------------------------
Analise post hoc:
------------------------------------------

Residuos ajustados standardizados corrigidos por momento (MCSTARs):
            Tabagista NaoTabagista
Etilista     6.963186    -6.963186
NaoEtilista -6.963186     6.963186

|MCSTAR critico| (alfaBonferroni=5%/1) = 1.959964

------------------------------------------
Analise de significancia pratica:
------------------------------------------

V de Cramer = 0.3319569 
Grau intermediario de dependencia entre as duas variaveis nominais

Teste de razao de chances (OR) robusto:

    Fisher's Exact Test for Count Data

data:  TC
p-value = 0.0006368
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
  1.694576 10.331999
sample estimates:
odds ratio 
  4.107342 

A novidade está na sessão Residuos ajustados standardizados corrigidos por momento (MCSTARs). Analise pelos valores positivos acima do valor assinalado como [MCSTAR critico]: são escores \(z\) relativos a quanto cada célula contribuiu para o qui-quadrado estimado. Neste exemplo, a diagonal principal é a responsável pela rejeição de \(H_0\), indicando associação concordante entre etilismo e tabagismo.


A única “sujeira” deste código é:

New names:
* `` -> …1

Isto desaparece se colocarmos algo, que será desprezado, na célula A1 da planilha. Experimente rodar o código com a planilha tabagismo_e_etilismo_2x2_v2.xlsx.

Note, também, o uso das linhas:

options(warn=-1)
[…]
options(warn=0)

que são as responsáveis por inibir os warnings do código (só as use depois que seu código estiver correto; warnings são importantes para evitar erros enquanto desenvolve um script).

Qui-quadrado mínimo

Segundo este autor, a correção de Yates não deve ser usada:



Note a relação entre esta forma de calcular o qui-quadrado e o odds-ratio:

Qui-quadrado para variáveis ordinais

O teste qui-quadrado tradicional trata os dois fatores com categorias nominais. Há muitos casos, porém, em que estes fatores são ordinais. Embora não seja incorreto utilizar o qui-quadrado (toda variável ordinal é também nominal), alguma informação está sendo perdida.

Existe uma medida, o Gama de Goodman-Kruskal, que considera ambas as variáveis como categóricas ordinais, potencialmente melhorando o resultado do teste. É uma estatistica nao-parametrica que mede a força de associação em dados cruzados em tabela, com as hipóteses:

\[H_o: \gamma = 0\] \[H_o: \gamma \ne 0\]

Neste estudo, Goodman (1979) desenvolve a análise com variáveis ordinais.

Os dados desta tabela são, originalmente, de:

Duncan, O.D. (1979), How Destination Depends on Origin in the Occupational Mobility Table. American Journal of Sociology, 84, 793-803.

As categorias são as seguintes:

    1. Professional and high administrative
    1. Managerial and executive
    1. Inspectional, supervisory, and other nonmanuar (high grade)
    1. Inspectional, supervisory, and other nonmanual (lower grade)
  1. (V)(a) Routine grades of nonmanual
  2. (V)(b) Skilled manual
    1. Semiskilled manual
    1. Unskilled manual

Portanto, a variável de status ocupacional é, claramente, qualitativa ordinal.


O repositório do R tem mais que pacotes. Há vários bancos de dados armazenados nele, para uso em exemplos. Aqui usaremos os dados deste estudo, disponíveis no pacote “dataset”.

Implementado em TCrXc_StatusOcupacional.R:

# TCrXc_StatusOcupacional.R

library(datasets)
library(DescTools)

eiras.show.MCSTAR <- function (MCSTAR, zcrit)
{
  MCSTARtxt <- round(MCSTAR,3)
  for (r in 1:nrow(MCSTAR))
  {
    for (c in 1:nrow(MCSTAR))
    {
      if (MCSTAR[r,c] > zcrit)
      {
        MCSTARtxt[r,c] <- paste("#",MCSTARtxt[r,c],sep="")
      } else
      {
        MCSTARtxt[r,c] <- paste(" ",MCSTARtxt[r,c],sep="")
      }
    }
  }
  return(MCSTARtxt)
}

# desabilita warnings
options(warn=-1)

cat("\n------------------------------------------\n")
cat("Dados:")
cat("\n------------------------------------------\n")
# Duncan, O.D. (1979), How Destination Depends on Origin in the 
# Occupational Mobility Table. American Journal of Sociology 84, 793-803.
tcrxc <- datasets::occupationalStatus
# para recolocar as categorias originais
categories <- c("I","II","III","IV","Va","Vb","VI","VII")
rownames(tcrxc) <- categories
colnames(tcrxc) <- categories
prmatrix(tcrxc)

cat("\n------------------------------------------\n")
cat("Analise de significancia estatistica:")
cat("\n------------------------------------------\n")
alfa <- 0.05

cat("\nTeste qui-quadrado de Pearson exato:\n")
print(res <- chisq.test(tcrxc,simulate.p.value=TRUE,B=1e6))
N <- sum(res$observed)
nL <- nrow(tcrxc) # ou dim(TC)[1]
nC <- ncol(tcrxc) # ou dim(TC)[2]
df <- (nL-1)*(nC-1)
cat("X^2 critico de 95% = ",qchisq(p=1-alfa,df=df), "\n", sep="")
X2 <- res$statistic # estatistica de teste qui-quadrado
cat("Graus de liberdade (nao fornecidos por bootstrapping): ", df, "\n", sep="")
cat("Heuristica de significancia (rej. H0 se X^2/gl > 2): ", X2/df, "\n", sep="")

cat("\n------------------------------------------\n")
cat("Analise post hoc:")
cat("\n------------------------------------------\n")
cat("\nResiduos ajustados standarizados corrigidos por momento (MCSTARs):\n")
STAR <- res$stdres # STandardized Adjusted Residual (STAR)
MCSTAR <- STAR/(sqrt((1-1/nL)*(1-1/nC))) # Moment-correct (MCSTAR)
alfaBonf <- alfa/df
zcrit <- abs(qnorm(alfaBonf/2))
print(eiras.show.MCSTAR(MCSTAR,zcrit))
cat("\n|MCSTAR critico| (alfaBonferroni=5%/",df,") = ",zcrit,"\n", sep="")
phi2 <- X2/N # phi ou w de Cohen
Dim <- min(nL,nC)-1 # dimensoes da TC: dim = max(phi^2)

cat("\n------------------------------------------\n")
cat("Analise de significancia pratica:")
cat("\n------------------------------------------\n")
V <- sqrt(phi2/Dim) # V ou C de Cramer
if (0 <= V & V < 0.1) {gV <- "minimo"}
if (0.1 <= V & V < 0.3) {gV <- "pequeno"}
if (0.3 <= V & V < 0.5) {gV <- "intermediario"}
if (0.5 <= V & V <= 1.0) {gV <- "grande"}
cat("\nV de Cramer =", V, "\nGrau", gV, 
    "de dependencia entre as duas variaveis nominais\n")

# Goodman, L. A. (1979) Simple Models for the Analysis of Association 
# in Cross-Classifications having Ordered Categories. 
# J. Am. Stat. Assoc., 74 (367), 537–552.
gama <- DescTools::GoodmanKruskalGamma(tcrxc)
cat("\nGama de Goodman-Kruskal = ", gama,"\n", sep="")

# EDWARDES, MD & BALTZAN, M (2000) The generalization of the odds ratio,
# rik ratio and risk difference to rxk tables.
# Statistics in Medicine 19:1901-14.
ORg <- (1+gama)/(1-gama)
cat("\nOR generalizado = (1+",gama,")/(1-",gama,") = ", ORg,"\n", sep="")

# reabilita warnings
options(warn=0)

A saída é:


------------------------------------------
Dados:
------------------------------------------
     I II III  IV Va  Vb  VI VII
I   50 19  26   8  7  11   6   2
II  16 40  34  18 11  20   8   3
III 12 35  65  66 35  88  23  21
IV  11 20  58 110 40 183  64  32
Va   2  8  12  23 25  46  28  12
Vb  12 28 102 162 90 554 230 177
VI   0  6  19  40 21 158 143  71
VII  0  3  14  32 15 126  91 106

------------------------------------------
Analise de significancia estatistica:
------------------------------------------

Teste qui-quadrado de Pearson exato:

    Pearson's Chi-squared test with simulated p-value (based on 1e+06
    replicates)

data:  tcrxc
X-squared = 1416, df = NA, p-value = 1e-06

X^2 critico de 95% = 66.33865
Graus de liberdade (nao fornecidos por bootstrapping): 49
Heuristica de significancia (rej. H0 se X^2/gl > 2): 28.89877

------------------------------------------
Analise post hoc:
------------------------------------------

Residuos ajustados standarizados corrigidos por momento (MCSTARs):
      destination
origin I       II      III     IV      Va      Vb      VI      VII    
   I   #28.022 #6.466  #4.851   -2.711  -0.804  -7.091  -4.336  -4.284
   II  #6.535  #15.194 #6.477   -0.475  0.201   -6.217  -4.43   -4.437
   III  0.706  #6.01   #7.195  #3.979   2.782   -3.966  -6.129  -4.134
   IV   -1.369  -0.926  1.7    #6.772   0.826   0.847   -3.453  -5.132
   Va   -1.436  0.409   -0.87   0.701  #5.188   -1.363  0.388   -1.982
   Vb   -6.546  -6.397  -3.505  -1.856  -0.703 #7.926   0.031   1.551 
   VI   -4.57   -4.075  -4.744  -3.41   -2.462  0.329  #9.978   2.718 
   VII  -4.152  -4.315  -4.744  -3.427  -2.901  -0.678 #4.169  #11.153

|MCSTAR critico| (alfaBonferroni=5%/49) = 3.284839

------------------------------------------
Analise de significancia pratica:
------------------------------------------

V de Cramer = 0.2404799 
Grau pequeno de dependencia entre as duas variaveis nominais

Gama de Goodman-Kruskal = 0.4209053

OR generalizado = (1+0.4209053)/(1-0.4209053) = 2.453667

É praticamente o mesmo código utilizado para o exemplo de tabagismo e etilismo, com algumas diferenças:

  • adicionamos uma função, eiras.show.MCSTAR() para exibir a matrix resultante de MCSTAR com menor número de casas decimais e símbolos adicionais (“#”) para que os valores acima de \(z\) crítico sejam mais facilmente localizados.
  • em uma tabela com 8 linhas e 8 colunas não existe definição para o odds ratio. No lugar deste, aparece o gama de Godman-Kruskal e um cálculo de OR generalizado.

Apenas para comparação, o teste qui-quadrado comum produziria o seguinte:

source("TCrXc_StatusOcupacional.R")

------------------------------------------
Dados:
------------------------------------------
     I II III  IV Va  Vb  VI VII
I   50 19  26   8  7  11   6   2
II  16 40  34  18 11  20   8   3
III 12 35  65  66 35  88  23  21
IV  11 20  58 110 40 183  64  32
Va   2  8  12  23 25  46  28  12
Vb  12 28 102 162 90 554 230 177
VI   0  6  19  40 21 158 143  71
VII  0  3  14  32 15 126  91 106

------------------------------------------
Analise de significancia estatistica:
------------------------------------------

Teste qui-quadrado de Pearson exato:

    Pearson's Chi-squared test with simulated p-value (based on 1e+06
    replicates)

data:  tcrxc
X-squared = 1416, df = NA, p-value = 1e-06

X^2 critico de 95% = 66.33865
Graus de liberdade (nao fornecidos por bootstrapping): 49
Heuristica de significancia (rej. H0 se X^2/gl > 2): 28.89877

------------------------------------------
Analise post hoc:
------------------------------------------

Residuos ajustados standarizados corrigidos por momento (MCSTARs):
      destination
origin I       II      III     IV      Va      Vb      VI      VII    
   I   #28.022 #6.466  #4.851   -2.711  -0.804  -7.091  -4.336  -4.284
   II  #6.535  #15.194 #6.477   -0.475  0.201   -6.217  -4.43   -4.437
   III  0.706  #6.01   #7.195  #3.979   2.782   -3.966  -6.129  -4.134
   IV   -1.369  -0.926  1.7    #6.772   0.826   0.847   -3.453  -5.132
   Va   -1.436  0.409   -0.87   0.701  #5.188   -1.363  0.388   -1.982
   Vb   -6.546  -6.397  -3.505  -1.856  -0.703 #7.926   0.031   1.551 
   VI   -4.57   -4.075  -4.744  -3.41   -2.462  0.329  #9.978   2.718 
   VII  -4.152  -4.315  -4.744  -3.427  -2.901  -0.678 #4.169  #11.153

|MCSTAR critico| (alfaBonferroni=5%/49) = 3.284839

------------------------------------------
Analise de significancia pratica:
------------------------------------------

V de Cramer = 0.2404799 
Grau pequeno de dependencia entre as duas variaveis nominais

Gama de Goodman-Kruskal = 0.4209053

OR generalizado = (1+0.4209053)/(1-0.4209053) = 2.453667

Especificamente para tabelas 2x2 existe o \(Q\) de Yule, dado por \[Q = {{ad-bc} \over {ad+bc}}\] Assim como o Gama de Goodman-Kruskal, varia de \(-1\) a \(+1\), mas como se aplica a tabelas 2x2, relaciona-se com o odds-ratio: \[Q = {{OR-1} \over {OR+1}}\] Uma função está implementada em psych::Yule(x, Y=FALSE), que devolve o \(Q\) de Yule quando \(x\) é uma tabela de frequências absolutas 2x2.

Qui-quadrado para tendência de Cochran-Armitage

Para situações em que uma variável é ordinal e a outra é nominal, por exemplo quando há um modelo de dose-resposta com um desfecho dicotômico associado a uma exposição que tem mais que dois níveis, é preferível utilizar o teste conhecido como Cochran-Armitage test for trend.

exemplo: espessura da dobra cutânea e menarca

Os pesquisadores Beckles et al. (1985) estudaram a associação entre a variável de a faixa etária da menarca (\(<\) 12 anos e $$12 anos) e o fator a espessura da dobra cutânea do tríceps (pequena, intermediária e grande).
Beckles et al. (1985) International Journal of Obesity 9:127-35, apud Kirkwood & Sterne (2003): Chapter 17: Chi-squared tests for 2x2 and larger contingency tables.

Caso ambas as variáveis sejam consideradas nominais, o teste adequado é o qui-quadrado de Pearson. Caso a variável considerada como exposição seja ordinal e o desfecho seja dicotômico (tabela \(K\)x\(2\)), o teste qui-quadrado de Cochran Armitage pode ser usado.

Para comparação, ambos os testes são executados em Teste qui-quadrado de Cochran-Armitage.R:

# Teste qui-quadrado de Cochran-Armitage.R

# Cochran-Armitage test for trend
# Menarca precoce & espessura da dobra cutânea do tríceps  
# Fonte:  Beckles et al. (1985) International Journal of Obesity 9:127-35,
# apud Kirkwood & Sterne (2003): Chapter 17: Chi-squared tests for 2x2 and 
# larger contingency tables. 
library(readxl)
library(DescTools)

# desabilita warnings
options(warn=-1)

# Tabela de Dancey&Reidy (2019, p. 275)
TCtmp <- read_excel("menarca.xlsx")
TC <- data.matrix(TCtmp)
rownames(TC) <-as.character(unlist(TCtmp[,1]))
TC <- TC[,-1]
remove(TCtmp)
cat("\n------------------------------------------\n")
cat("Dados:")
cat("\n------------------------------------------\n")
print(TC)

cat("\n------------------------------------------\n")
cat("Analise de significancia estatistica:")
cat("\n------------------------------------------\n")
alfa <- 0.05

cat("\nTeste qui-quadrado de Pearson exato:\n")
print(res <- chisq.test(TC,simulate.p.value=TRUE,B=1e6))
cat("X^2 critico de 95% = ",qchisq(p=1-alfa,df=1), "\n", sep="")
N <- sum(res$observed)
nL <- nrow(TC) # ou dim(TC)[1]
nC <- ncol(TC) # ou dim(TC)[2]
df <- (nL-1)*(nC-1)
X2 <- res$statistic # estatistica de teste qui-quadrado
cat("Graus de liberdade (nao fornecidos por bootstrapping): ", df, "\n", sep="")
cat("Heuristica de significancia (rej. H0 se X^2/gl > 2): ", X2/df, "\n", sep="")

cat("\nTeste qui-quadrado de Cochran-Armitage:\n")
res <- DescTools::CochranArmitageTest(TC)
print(res)
x2ca <- res$statistic^2
p <- res$p.value
cat("\nX-squared trend =",x2ca, ", df = 1, ", "p-value =", p,"\n")

cat("\n------------------------------------------\n")
cat("Analise de significancia pratica:")
cat("\n------------------------------------------\n")
phi2 <- X2/N # phi ou w de Cohen
Dim <- min(nL,nC)-1 # dimensoes da TC: dim = max(phi^2)
V <- sqrt(phi2/Dim) # V ou C de Cramer
if (0 <= V & V < 0.1) {gV <- "minimo"}
if (0.1 <= V & V < 0.3) {gV <- "pequeno"}
if (0.3 <= V & V < 0.5) {gV <- "intermediario"}
if (0.5 <= V & V <= 1.0) {gV <- "grande"}
cat("\nV de Cramer =", V, "\nGrau", gV, 
    "de dependencia entre as duas variaveis nominais\n")

# reabilita warnings
options(warn=0)

que produz:


------------------------------------------
Dados:
------------------------------------------
              menor12 maiorigual12
Pequena            15          156
Intermediaria      29          197
Grande             36          186

------------------------------------------
Analise de significancia estatistica:
------------------------------------------

Teste qui-quadrado de Pearson exato:

    Pearson's Chi-squared test with simulated p-value (based on 1e+06
    replicates)

data:  TC
X-squared = 4.7594, df = NA, p-value = 0.09404

X^2 critico de 95% = 3.841459
Graus de liberdade (nao fornecidos por bootstrapping): 2
Heuristica de significancia (rej. H0 se X^2/gl > 2): 2.379692

Teste qui-quadrado de Cochran-Armitage:

    Cochran-Armitage test for trend

data:  TC
Z = 2.1783, dim = 3, p-value = 0.02938
alternative hypothesis: two.sided


X-squared trend = 4.744926 , df = 1,  p-value = 0.02938481 

------------------------------------------
Analise de significancia pratica:
------------------------------------------

V de Cramer = 0.08768596 
Grau minimo de dependencia entre as duas variaveis nominais

Compare os valores \(p\) de ambos os testes. Cochran-Armitage, ao considerar a espessura da prega cutânea como “dose”, ordinal, indica significância estatística onde o teste convencional não detectava.

Cochran-Armitage generalizado

Para desfecho com três ou mais categorias (os desfechos são tratados como variáveis nominais) existe função no pacote coin::chisq.test

Um exemplo está em Teste qui-quadrado de Cochran-Armitage generalizado.R:

library(coin)
library(rcompanion)
# The chisq_test function in the coin package can be used to conduct a test of 
# association for a contingency table with one ordered nominal variable and 
# one non-ordered nominal variable.  
# The Cochran–Armitage test is a special case of this where the 
# non-ordered variable has only two categories.
# The scores option is used to indicate which variable should be treated as 
# ordered, and the spacing of the levels of this variable.
Job <- matrix(c(1,2,1,0, 3,3,6,1, 10,10,14,9, 6,7,12,11), 4, 4,
              dimnames = list(income = c("< 15k", "15-25k", "25-40k", "> 40k"),
                              satisfaction = c("VeryD", "LittleD", 
                                               "ModerateS", "VeryS")))
Job <- as.table(Job)
print(coin::chisq_test(Job,  scores = list("income"= c(1, 2, 3, 4))))
# Teste post hoc
print(ph <- rcompanion::pairwiseOrdinalIndependence(Job,compare="column"))
try(rcompanion::cldList(p.value ~ Comparison, data = ph, threshold = 0.05))

que produz:

Loading required package: survival

Attaching package: 'survival'
The following object is masked from 'package:epitools':

    ratetable

    Asymptotic Generalized Pearson Chi-Squared Test

data:  satisfaction by
     income (< 15k < 15-25k < 25-40k < > 40k)
chi-squared = 3.1362, df = 3, p-value = 0.3711

           Comparison p.value p.adjust
1     VeryD : LittleD   0.451    0.541
2   VeryD : ModerateS   0.351    0.526
3       VeryD : VeryS   0.161    0.526
4 LittleD : ModerateS   0.698    0.698
5     LittleD : VeryS   0.242    0.526
6   ModerateS : VeryS   0.271    0.526

Qui-quadrado de Mantel-Haenszel

Este teste utiliza uma tabela tridimensional. Neste exemplo temos a contagem de mortes entre fumantes e não fumantes estratificadas por faixa etária. Esta versão de Qui-quadrado faz o teste global e, depois, estratifica pela terceira variável (faixa etária, neste exemplo).

Faixa Etária 18-24 25-34 35-44 45-54 55-64 65-74 75
Tabagista Desfecho
Sim Morta 2 3 14 27 51 29 13
Viva 53 121 95 103 64 7 0
Não Morta 1 5 7 12 40 101 64
Viva 61 152 114 66 81 28 0

Os dados precisam ser reorganizados para que construamos uma tabela tridimensional requerida pela funções relacionadas a este teste. Veja como os dados foram acomodados na planila fumo_e_faixaetaria.xlsx.

Implementamos o teste em Teste Qui-quadrado de Mantel-Haenszel.R:

# Teste Qui-quadrado de Mantel-Haenszel.R

# desabilita warnings
options(warn=-1)

# Tabela de contingencia 2x2 segmentada  
# Teste robusto da razao de chances (OR) de Mantel-Haenszel

# Relacao entre tabagismo e sobrevivência em 20 anos (1974-1994) 
# em 1.134 mulheres adultas do Reino Unido
# Delineamento: coorte
# Fonte:  APPLETON, D. R. et al. (1996) Ignoring a covariate:
# An example of Simpson's paradox. The American Statistician,
# 50(4): 340-1.

# A aspa inicial TEM que comecar na primeira coluna da linha
# e os espacamentos distintos dos dessa tabela podem causar 
# problemas de na geracao da tabela horizontalizada (ftable).
library(vcd)
library(rcompanion)
library(coin)
library(readxl)
library(tidyverse)
TCtmp <- read_excel("fumo_e_faixaetaria.xlsx")
TCtmp <- gather(TCtmp, var, Contagem, -Idade)
TCtmp <- separate(TCtmp, var, c('Exposicao','Desfecho')) 
TCS <- TCtmp %>% xtabs(Contagem ~ Exposicao + Desfecho + Idade, .)
print(TCS)

# Mantel-Haenszel test of the null that two nominal variables 
# are conditionally independent in each stratum, 
# assuming that there is no three-way interaction.
# Woolf test for homogeneity of ORs across strata on 2 x 2 x k tables:
# If significant, M-H test is not appropriate.
print(ftable(TCS)) # Display a flattened table (tabela horizontalizada)

cat("\nHomogeneidade dos ORs:\n")
print(vcd::woolf_test(TCS)) 

# Teste para tabela 2x2xK
cat("\nQui-quadrado de Mantel-Haenszel:\n")
print(mantelhaen.test(TCS, exact=TRUE)) # Teste exato de OR de Mantel-Haenszel

# Teste para tabela LxCxK
cat("\nTeste de Cochran-Mantel-Haenszel (coin::cmh_test):\n")
print(coin::cmh_test(TCS, distribution = approximate(nresample = 1e6)))

cat("\nAnalise estratificada por idade:\n")
print(rcompanion::groupwiseCMH(TCS,
             group   = 3,
             fisher  = TRUE,
             gtest   = FALSE,
             chisq   = FALSE,
             method  = "fdr",
             correct = "none",
             digits  = 3))
cat("\n")
print(ors <- vcd::oddsratio(TCS, log=FALSE)) # Show OR for each 2x2
cat("\n")
lnors <- vcd::oddsratio(TCS, log=TRUE) # Show ln(OR) for each 2x2
print(lnors)
print(summary(lnors))
print(confint(lnors))
plot(lnors, xlab = "Faixa Etaria", ylab = "ln(OR)")

# reabilita warnings
options(warn=0)

que gera a seguinte saída:

Loading required package: grid

Attaching package: 'vcd'
The following object is masked from 'package:epitools':

    oddsratio
── Attaching packages ──────────────────────────────────────────────────── tidyverse 1.3.0 ──
✓ ggplot2 3.3.0     ✓ purrr   0.3.3
✓ tibble  3.0.1     ✓ dplyr   0.8.5
✓ tidyr   1.0.2     ✓ stringr 1.4.0
✓ readr   1.3.1     ✓ forcats 0.5.0
── Conflicts ─────────────────────────────────────────────────────── tidyverse_conflicts() ──
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()
, , Idade = 18-24

              Desfecho
Exposicao      Morta Viva
  NaoTabagista     1   61
  Tabagista        2   53

, , Idade = 25-34

              Desfecho
Exposicao      Morta Viva
  NaoTabagista     5  152
  Tabagista        3  121

, , Idade = 35-44

              Desfecho
Exposicao      Morta Viva
  NaoTabagista     7  114
  Tabagista       14   95

, , Idade = 45-54

              Desfecho
Exposicao      Morta Viva
  NaoTabagista    12   66
  Tabagista       27  103

, , Idade = 55-64

              Desfecho
Exposicao      Morta Viva
  NaoTabagista    40   81
  Tabagista       51   64

, , Idade = 65-74

              Desfecho
Exposicao      Morta Viva
  NaoTabagista   101   28
  Tabagista       29    7

, , Idade = 75+

              Desfecho
Exposicao      Morta Viva
  NaoTabagista    64    0
  Tabagista       13    0

                      Idade 18-24 25-34 35-44 45-54 55-64 65-74 75+
Exposicao    Desfecho                                              
NaoTabagista Morta              1     5     7    12    40   101  64
             Viva              61   152   114    66    81    28   0
Tabagista    Morta              2     3    14    27    51    29  13
             Viva              53   121    95   103    64     7   0

Homogeneidade dos ORs:

    Woolf-test on Homogeneity of Odds Ratios (no 3-Way assoc.)

data:  TCS
X-squared = 3.2061, df = 6, p-value = 0.7826


Qui-quadrado de Mantel-Haenszel:

    Exact conditional test of independence in 2 x 2 x k tables

data:  TCS
S = 230, p-value = 0.01591
alternative hypothesis: true common odds ratio is not equal to 1
95 percent confidence interval:
 0.4538410 0.9355508
sample estimates:
common odds ratio 
        0.6534856 


Teste de Cochran-Mantel-Haenszel (coin::cmh_test):

    Approximative Generalized Cochran-Mantel-Haenszel Test

data:  Desfecho by
     Exposicao (NaoTabagista, Tabagista) 
     stratified by Idade
chi-squared = 5.8443, p-value = 0.01584


Analise estratificada por idade:

  Group   Test p.value adj.p
1 18-24 Fisher  0.6000 1.000
2 25-34 Fisher  1.0000 1.000
3 35-44 Fisher  0.0706 0.291
4 45-54 Fisher  0.3650 0.852
5 55-64 Fisher  0.0830 0.291
6 65-74 Fisher  1.0000 1.000
7   75+ Fisher  1.0000 1.000

 odds ratios for Exposicao and Desfecho by Idade 

    18-24     25-34     35-44     45-54     55-64     65-74       75+ 
0.5219512 1.2519906 0.4314109 0.7074504 0.6223718 0.9054416 4.7777778 

log odds ratios for Exposicao and Desfecho by Idade 

      18-24       25-34       35-44       45-54       55-64       65-74 
-0.65018114  0.22473479 -0.84069420 -0.34608770 -0.47421763 -0.09933253 
        75+ 
 1.56397554 

z test of coefficients:

       Estimate Std. Error z value Pr(>|z|)  
18-24 -0.650181   1.049580 -0.6195  0.53561  
25-34  0.224735   0.694493  0.3236  0.74624  
35-44 -0.840694   0.470642 -1.7863  0.07406 .
45-54 -0.346088   0.375584 -0.9215  0.35681  
55-64 -0.474218   0.268109 -1.7687  0.07694 .
65-74 -0.099333   0.460621 -0.2156  0.82926  
75+    1.563976   2.022270  0.7734  0.43930  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

           2.5 %     97.5 %
18-24 -2.7073204 1.40695808
25-34 -1.1364462 1.58591573
35-44 -1.7631351 0.08174672
45-54 -1.0822181 0.39004270
55-64 -0.9997024 0.05126713
65-74 -1.0021328 0.80346776
75+   -2.3996018 5.52755287